#! /usr/bin/env python
#
"""Classes and functions to deal with fatgraphs.
"""
__docformat__ = 'reStructuredText'

from native cimport *

## logging subsystem

import logging

## stdlib imports

import sys
import operator
import itertools
import types
import weakref


## application-local imports

from cache cimport Cacheable, cache_id
from combinatorics cimport InplacePermutationIterator, SetProductIterator, Permutation, PermutationIterator
from cycliclist cimport CyclicList
from iterators cimport BufferingIterator, IndexedIterator, Iterator
from loadsave cimport load,save
from runtime import runtime
from utils cimport lconcat, sign, ltranslate


## aux functions (Cython doesn't support nested `def`s)

# ----- Fatgraph.isomorphisms()
cdef inline list starting_vertices(Fatgraph graph):
    """
    Return the pair `(valence, vertices)`, which minimizes the product
    of valence with the number of vertices of that valence.  Second
    field `vertices` is really the *indices* of the vertices in
    `graph`.

    Examples::

      >>> g = Fatgraph([Vertex([0,1,2]), Vertex([0,2,3]), Vertex([1,6,6]), Vertex([3,4,4,5,5,])])
      >>> starting_vertices(g)
      (5, [3])
    """
    cdef int val
    cdef int n
    cdef int n_
    val = (max(graph.vertex_valences()))
    vs = None
    n = (len(graph.vertices)+1)
    for (val_, vs_) in graph.valence_spectrum().iteritems():
        n_ = (len(vs_))
        if (n_*val_ < n*val) \
               or (n_*val_ == n*val and val_<val):
            val = val_
            vs = vs_
            n = n_
    return list([val, vs])

cdef inline bint compatible(Vertex v1, Vertex v2):
    """Return `True` if vertices `v1` and `v2` are compatible.
    (i.e., same valence and number of loops - one *could* be
    mapped onto the other.)
    """
    if len(v1) == len(v2) and v1.num_loops == v2.num_loops:
        return True
    else:
        return False

cdef inline list admissible_vertex_mappings(Vertex v, Fatgraph g, list ixs):
    """Iterate over all (indices of) vertices in `g`, which
    `v` *could* be mapped to (that is, the destination vertex
    matches `v` in valence and number of loops.

    Third argument `ixs` restricts the search to the given
    subset of vertex indices in `g`.
    """
    return [i for i in ixs if compatible(v, g.vertices[i])]

cdef class CannotExtendMap(BaseException):
    """Exception raised by `extend_map` on failure to extend a
    partial map.
    """
    pass

cdef inline extend_map(m, Fatgraph g1, int i1, int r, Fatgraph g2, int i2):
    """Extend map `m` by mapping the `i1`-th vertex in `g1` to
    the `i2`-th vertex in `g2` (and rotating the source vertex
    by `r` places leftwards).  Return the extended map `(pv,
    rot, pe)`.

    The partial map `m` is represented as a triple `(pv, rot,
    pe)` as in `Fatgraph.isomorphism` (which see), with the
    additional proviso that unassigned items in `rot` are
    represented by `None`.
    """
    (pv, rots, pe) = m
    v1 = g1.vertices[i1]
    v2 = g2.vertices[i2]
    if not compatible(v1, v2):
        raise CannotExtendMap

    # XXX: rotation has to be >=0 for the [r:r+..] shift below to work
    if r < 0:
        r += len(v2)

    if pv.has_key(i1):
        if pv[i1] != i2 or (rots[i1] - r) % len(v2) != 0:
            raise CannotExtendMap
        else:
            # this pair has already been added
            return m

    pv[i1] = i2
    rots[i1] = r

    # rotating `v1` leftwards is equivalent to rotating `v2` rightwards...
    v2 = v2[r:r+len(v2)]
    if not pe.extend(v1, v2):
        raise CannotExtendMap

    return (pv, rots, pe)

cdef inline list neighbors(m, Fatgraph g1, int i1, Fatgraph g2, int i2):
    """List of vertex-to-vertex mappings that extend map
    `m` in the neighborhood of of `i1` (in the domain) and
    `i2` (in the codomain).

    Return a list of triplets `(src, dst, rot)`, where:
       * `src` is the index of a vertex in `g1`,
         connected to `i1` by an edge `x`;
       * `dst` is the index of a vertex in `g2`,
         connected to `i2` by the image (according to `m`)
         of edge `x`;
       * `rot` is the rotation to be applied to `g1.vertices[src]`
         so that edge `x` and its image appear
         at the same index position;

    XXX: prune vertices that are already mapped by `m`?
    """
    cdef int dst_i
    cdef int dst_v
    cdef int src_i
    cdef int src_v
    cdef int x
    cdef list dst_endpoints_i
    cdef list dst_endpoints_v
    cdef list result
    cdef list src_endpoints_i
    cdef list src_endpoints_v

    result = list()
    for x in g1.vertices[i1]:
            src_endpoints_v = g1.endpoints_v[x]
            # ignore loops
            if src_endpoints_v[0] == src_endpoints_v[1]:
                continue # to next `x`
            src_endpoints_i = g1.endpoints_i[x]
            src_v = src_endpoints_v[0] if (src_endpoints_v[1] == i1) else src_endpoints_v[1]
            # ignore vertices that are already in the domain of `m`
            if src_v in m[0]:
                continue # to next `x`
            src_i = src_endpoints_i[0] if (src_endpoints_v[1] == i1) else src_endpoints_i[1]
            dst_endpoints_v = g2.endpoints_v[m[2][x]]
            dst_endpoints_i = g2.endpoints_i[m[2][x]]
            dst_v = dst_endpoints_v[0] if (dst_endpoints_v[1] == i2) else dst_endpoints_v[1]
            dst_i = dst_endpoints_i[0] if (dst_endpoints_v[1] == i2) else dst_endpoints_i[1]
            # array of (source vertex index, dest vertex index, rotation)
            result.append((src_v, dst_v, dst_i-src_i))
    return result

# ------- end of Fatgraph.isomorphisms()


# ------- Fatgraph.compute_boundary_cycles()

cdef inline first_unused_corner(list vertices, list corners):
    cdef int v
    cdef int i
    cdef corner
    # fast-forward to the first unused corner
    for v in xrange(len(vertices)):
        for i in xrange(len(vertices[v])):
            corner = corners[v][i]
            if corner is not None:
                return (v, i, corner)
    return (v, i, None)

# ------- end of Fatgraph.compute_boundary_cycles()



## main

cdef class BoundaryCycle(frozenset):
    """A boundary cycle of a Fatgraph.

    Boundary cycles are a cyclic sequence of 'corners': a corner
    consists of a vertex `v` and (an unordered pair of) two
    consecutive indices (in the cyclic order at `v`, so, either `j
    == i+1` or `i` and `j` are the starting and ending indices).

    Two boundary cycles are equal if they comprise the same
    corners.
    """

    cpdef BoundaryCycle set_graph(self, Fatgraph graph):
        # use a weakref so not to create a reference cycle and ease GC
        self.graph = weakref.proxy(graph)
        return self
        

    cpdef contract_bcy(self, vi1, vi2, Fatgraph graph):
        """Return a new `BoundaryCycle` instance, image of the
        topological map that contracts the edge with endpoints
        `(v1,i1)` and `(v2,i2)` that are passed as first and
        second argument.

        Optional third argument `graph` is passed unchanged to the
        `BoundaryCycle` constructor.

        XXX: return `self` if neither `v1` nor `v2` are contained
        here.
        """
        cdef int i1
        cdef int i2
        cdef int l1
        cdef int l2
        cdef int pos1
        cdef int pos2
        cdef int v1
        cdef int v2
        cdef list new_bcy
        cdef object corner # tuple
        
        (v1, pos1) = vi1
        (v2, pos2) = vi2
        l1 = len(self.graph.vertices[v1])
        l2 = len(self.graph.vertices[v2])
        new_bcy = list()
        for corner in self:
            if corner[0] == v1:
                if pos1 == corner[1]:
                    # skip this corner, keep only one of the
                    # corners limited by the contracted edge
                    continue
                else: 
                    i1 = (corner[1] - pos1 - 1) % l1
                    i2 = (corner[2] - pos1 - 1) % l1
                    assert (i1+1-i2) % l1 == 0 # i1,i2 denote successive indices
                    assert i1 != l1-1 # would collide with contracted corners from `v2`
                    new_bcy.append((v1, i1, i2))
            elif corner[0] == v2:
                if pos2 == corner[1]:
                    # skip this corner, keep only one of the
                    # corners limited by the contracted edge
                    continue
                if pos2 == corner[2]:
                    new_bcy.append((v1, l1+l2-3, 0))
                else:
                    i1 = l1-1 + ((corner[1] - pos2 - 1) % l2)
                    i2 = l1-1 + ((corner[2] - pos2 - 1) % l2)
                    assert (i1+1-i2) % l1 == 0 # i1,i2 denote successive indices
                    new_bcy.append((v1, i1, i2))
            elif corner[0] > v2:
                # shift vertices after `v2` one position down
                new_bcy.append((corner[0]-1, corner[1], corner[2]))
            else:
                # pass corner unchanged
                new_bcy.append(corner)
        return BoundaryCycle(new_bcy).set_graph(graph)

    cpdef transform(self, iso, Fatgraph graph):
        """Return a new `BoundaryCycle` instance, obtained by
        transforming each corner according to a graph isomorphism.
        """
        assert self.graph is not None
        (pv, rots, pe) = iso
        triples = []
        cdef int v
        cdef int i
        cdef int j
        cdef int l
        cdef int v_
        cdef int i_
        cdef int j_
        for (v, i, j) in self:
            l = (len(self.graph.vertices[v]))
            # create transformed triple 
            v_ = (pv[v])
            i_ = ((i + rots[v]) % l)
            j_ = ((j + rots[v]) % l)
            # ensure the contract is honored, that `j` is the
            # index _following_ `i` in the cyclic order
            if i_ == 0 and j_ == l:
                i_, j_ = j_, i_
            triples.append((v_, i_, j_))
        return BoundaryCycle(triples).set_graph(graph)



cdef class Vertex(CyclicList):
    """A (representative of) a vertex of a ribbon graph.

    A vertex is represented by the cyclically ordered list of its
    (decorated) edges.  The edge colorings may be accessed through a
    (read-only) sequence interface.

    At init time, the number of loops attached to this vertex is
    computed and stored in the `.num_loops` attribute::

        >>> Vertex([0,1,2]).num_loops
        0
        >>> Vertex([1,1,0,0]).num_loops
        2
    """

    def __init__(self, seq):
        CyclicList.__init__(self, seq)
        self.num_loops = len(self) - len(set(self))

    def __str__(self):
        return repr(self)



cdef class EqualIfIsomorphic(Cacheable):
    """Instances of this class will compare equal if there is an
    isomorphism mapping one to the other.
    """

    def __init__(self, invariants):
        self.invariants = invariants
        self._cache_eq = dict()
        # set this instance's _persistent_id
        Cacheable.__init__(self)


    def __cmp__(self, other):
        """Return 0 if `self` and `other` are isomorphic; -1 otherwise.
        """

        assert isinstance(other, type(self)), \
               "EqualIfIsomorphic.__eq__:" \
               " called with incompatible type arguments: `%s` and` %s`" % (type(self), type(other))

        # shortcuts
        if self is other:
            return 0
        if self.invariants != other.invariants:
            return -1

        # go the long way: try to find an explicit isomorphims
        # between `self` and `other`
        key = frozenset([self.cache_id(), other.cache_id()])
        if not self._cache_eq.has_key(key):
            self._cache_eq[key] = 0 if len(self.isomorphisms(other)) > 0 else -1
        return self._cache_eq[key]


    def __richcmp__(self, other, int op):
        assert type(self) == type(other), \
               "EqualIfIsomorphic.__richcmp__: " \
               " Cannot compare instances of different classes: %s with %s" \
               % (type(self), type(other))
        cdef int cmp_
        cmp_ = self.__cmp__(other)
        if op == 2: # ==
            return cmp_ == 0
        elif op == 3: # !=
            return cmp_ != 0
        else:
            raise NotImplementedError("EqualIfIsomorphic.__richcmp__: Can only test for equality.")



#: memoize isomorphisms between `Fatgraph` instances
cdef dict _cache_isomorphisms = dict()


cdef class Fatgraph(EqualIfIsomorphic):
    """A fully-decorated ribbon graph.

    Exports a (read-only) sequence interface, through which vertices
    can be accessed.

    Several attributes of this object are computed at init time and
    store characteristics of the fatgraph:

        `.boundary_cycles`
            List of the boundary components of this `Fatgraph` object.
            Each boundary component is represented by the list of
            (colored) edges.

        `.genus`
            The genus of the Riemann surface this fatgraph lies on.

        `.num_boundary_cycles`
            Number of boundary cycles of this `Fatgraph`.

        `.num_edges`
            The number of edges of this `Fatgraph` object.

        `.num_vertices`
            Number of vertices of this `Fatgraph` object.

        `.vertices`
            List of vertices of this `Fatgraph`; each one is an instance
            of the `Vertex` class; see `Fatgraph.__init__` for examples.

    Examples::
        >>> Fatgraph([Vertex([2,1,0]), Vertex([2,1,0])]).num_boundary_cycles
        1
        >>> Fatgraph([Vertex([2,1,0]), Vertex([2,0,1])]).num_boundary_cycles
        3

    Two `Fatgraph`s compare equal if they are isomorphic::

        >>> Fatgraph([Vertex([1,0,0,1])]) == Fatgraph([Vertex([1,1,0,0])])
        True

        >>> Fatgraph([Vertex([2,0,0]), Vertex([2,1,1])]) \
              == Fatgraph([Vertex([2,2,0]), Vertex([1,1,0])])
        True

        >>> Fatgraph([Vertex([2,0,1]), Vertex([2,0,1])]) \
              == Fatgraph([Vertex([2,1,0]), Vertex([2,0,1])])
        False

        >>> Fatgraph([Vertex([2,0,1]), Vertex([2,0,1])]) \
              == Fatgraph([Vertex([2,0,0]), Vertex([2,1,1])])
        False

        >>> Fatgraph([Vertex([2,0,0]), Vertex([2,1,1])]) \
              == Fatgraph([Vertex([1,1,0,0])])
        False
    """

    def __init__(self, g_or_vs, **kwargs):
        """Construct a `Fatgraph` instance, taking list of vertices.

        Argument `g_or_vs` can be either:

            - a sequence of `Vertex` class instances::  

                >>> g1 = Fatgraph([Vertex([2,0,1]), Vertex([2,1,0])])

            Note that the list of vertices is assigned, *not copied*
            into the instance variable.

        or:

          - another `Fatgraph` instance, in which case this acts as a
            copy-constructor (and any remaining arguments are
            ignored), returning a new `Fatgraph` instance::

              >>> g2 = Fatgraph(g1)
              >>> g2 is g1
              False
              >>> g2 == g1
              True

            The returned `Fatgraph` instance *shares* all attributes
            with the instance given as argument::

              >>> g2.vertices is g1.vertices
              True
              >>> g2.edge_numbering is g1.edge_numbering
              True
              >>> g2.endpoints_v is g1.endpoints_v
              True
              >>> g2.endpoints_i is g1.endpoints_i
              True

        """
        cdef int current_vertex_index
        cdef int edge_index_in_vertex
        cdef int edge
                
        # dispatch based on type of arguments passed
        if isinstance(g_or_vs, Fatgraph):
            # copy-constructor
            self.boundary_cycles = g_or_vs.boundary_cycles
            self.edge_numbering = g_or_vs.edge_numbering
            self.endpoints_i = g_or_vs.endpoints_i
            self.endpoints_v = g_or_vs.endpoints_v
            self.genus = g_or_vs.genus
            self.num_boundary_cycles = g_or_vs.num_boundary_cycles
            self.num_edges = g_or_vs.num_edges
            self.num_external_edges = g_or_vs.num_external_edges
            self.num_vertices = g_or_vs.num_vertices
            self.vertices = g_or_vs.vertices

        else: # initialize new instance

            #: list of vertices
            self.vertices = g_or_vs

            #: Number of edge colors
            self.num_edges = kwargs.get('num_edges',
                                        sum([len(v) for v in self.vertices]) / 2)

            #: Number of external (loose-end) edges
            self.num_external_edges = kwargs.get('num_external_edges', 0)

            #: Number of vertices  XXX: why is this settable with kwarg???
            self.num_vertices = kwargs.get('num_vertices', len(g_or_vs))

            if 'endpoints' in kwargs:
                (self.endpoints_v, self.endpoints_i) = kwargs.get('endpoints')
            else:
                #: Adjacency list of this graph.  For each edge, store a pair
                #: `(v1, v2)` where `v1` and `v2` are (indices of)
                #: endpoint vertices of an edge, and a corresponding pair
                #: `(i1, i2)` where `i1` and `i2` are indices of the given
                #: edge in vertices `v1` and `v2`.
                #: Each pair `(v, i)` represents a flag by the endpoint
                #: vertex and the index of the edge in the vertex.  (The
                #: vertex index alone is not enough for representing the
                #: edge arrow for loops.)
                self.endpoints_v = [ list() for dummy in xrange(self.num_edges) ]
                self.endpoints_i = [ list() for dummy in xrange(self.num_edges) ]
                for current_vertex_index in xrange(self.num_vertices):
                    for (edge_index_in_vertex, edge) \
                            in enumerate(self.vertices[current_vertex_index]):
                        assert edge in range(self.num_edges), \
                                   "Fatgraph.__init__:"\
                                   " edge number %d not in range 0..%d" \
                                   % (edge, self.num_edges)
                        self.endpoints_v[edge].append(current_vertex_index)
                        self.endpoints_i[edge].append(edge_index_in_vertex)

            ## Orientation is given by an ordering of the edges,
            ## which directly translates into an orientation of the
            ## associated cell.  
            if 'orientation' in kwargs:
                self.edge_numbering = kwargs.get('orientation')
            else:
                self.edge_numbering = [ x for x in xrange(self.num_edges) ]

            self.boundary_cycles = self.compute_boundary_cycles()
            self.num_boundary_cycles = len(self.boundary_cycles)

            # by Euler, K-L+n=2-2*g
            self.genus = (self.num_edges - self.num_vertices
                          - self.num_boundary_cycles + 2) / 2
        
        # caching
        self._cache_bcys = None
        self._cache_edge_orbits = None
        self._cache_edge_pair_orbits = None
        self._cache_isomorphisms = dict()
        self._cache_starting_vertices = None
        self._cache_valence_spectrum = None

        # before computing invariants, check that internal data
        # structures are in a consistent state
        assert self._ok()

        # used for isomorphism testing
        EqualIfIsomorphic.__init__(self, (
            self.num_vertices,
            self.num_edges,
            self.num_external_edges if self.num_external_edges > 0 else self.num_boundary_cycles,
            #self.vertex_valences(),
            #self.vertex_valence_distribution(),
            ))


    def _ok(self):
        """Perform coherency checks on internal state variables of
        `Fatgraph` instance and return `True` if they all pass.
        """
        assert self.num_edges > 0, \
               "Fatgraph `%s` has 0 edges." % (self)
        # check regular edges endpoints
        for (edge, ep_v, ep_i) in itertools.izip(itertools.count(),
                                       self.endpoints_v[:self.num_edges],
                                       self.endpoints_i[:self.num_edges]):
            assert isinstance(ep_v, list)  # Fatgraph.contract() updates this destructively
            assert isinstance(ep_i, list)  # Fatgraph.contract() updates this destructively
            assert len(ep_v) == 2
            assert len(ep_i) == 2
            assert isinstance(ep_v[0], int)
            assert isinstance(ep_v[1], int)
            assert isinstance(ep_i[0], int)
            assert isinstance(ep_i[1], int)
            assert (0 <= ep_v[0] < self.num_vertices)
            assert (0 <= ep_v[1] < self.num_vertices)
            assert (0 <= ep_i[0] < len(self.vertices[ep_v[0]]))
            assert (0 <= ep_i[1] < len(self.vertices[ep_v[1]])), \
                   "Fatgraph `%s`:"\
                   " invalid attachment indices `%s`" \
                   " for endpoints %s of regular edge %d"\
                   % (self, ep_i, ep_v, edge)
##                    "Fatgraph `%s` has invalid regular endpoints array `%s/%s`" \
##                    " invalid endpoints pair %s/%s for edge %d" \
##                    % (self, self.endpoints_v, self.endpoints_i,
##                       ep_v, ep_i, edge)
            assert (edge in self.vertices[ep_v[0]]), \
                    "Invalid endpoints %s for edge %d of graph `%s`" \
                    % (ep_v, edge, self)
            assert (edge in self.vertices[ep_v[1]]), \
                    "Invalid endpoints %s for edge %d of graph `%s`" \
                    % (ep_v, edge, self)
        # check external edges endpoints
        for (edge, ep_v, ep_i) in itertools.izip(itertools.count(),
                                       self.endpoints_v[self.num_edges:],
                                       self.endpoints_i[self.num_edges:]):
            assert isinstance(ep_v, list)  # Fatgraph.contract() updates this destructively
            assert isinstance(ep_i, list)  # Fatgraph.contract() updates this destructively
            xedge = -self.num_external_edges + edge
            assert (ep_v[1] is None)
            assert isinstance(ep_v[0], int)
            assert (ep_i[1] is None)
            assert isinstance(ep_i[0], int)
            assert (0 <= ep_v[0] < self.num_vertices)
            assert (0 <= ep_i[0] < len(self.vertices[ep_v[0]]))
##                    "Fatgraph `%s` has invalid external endpoints array: `%s/%s`" \
##                    % (self, self.endpoints_v, self.endpoints_i)
        # check that each edge occurs exactly two times in vertices
        cnt = [ 0 for x in xrange(self.num_edges + self.num_external_edges) ]
        for v in self.vertices:
            for edge in v:
                cnt[edge] += 1
        for edge, cnt in enumerate(cnt):
            if edge < self.num_edges:
                assert cnt == 2, \
                       "Regular edge %d appears in %d vertices" \
                       % (edge, cnt)
            else:
                assert cnt == 1, \
                       "External edge %d appears in %d vertices" \
                       % (edge, cnt)

        assert self.edge_numbering is not None
        
        return True


    def __repr__(self):
        if hasattr(self, 'num_external_edges') and self.num_external_edges > 0:
            return "Fatgraph(%s, num_external_edges=%d)" \
                   % (repr(self.vertices), self.num_external_edges)
        elif hasattr(self, 'vertices'):
            return "Fatgraph(%s)" % repr(self.vertices)
        else:
            # avoid PyDB crashing when stepping into `Fatgraph.__init__` code
            return "Fatgraph(<Initializing...>)"

    
    def __str__(self):
        return repr(self)


    cpdef list automorphisms(self):
        """Enumerate automorphisms of this `Fatgraph` object.

        See `.isomorphisms()` for details of how a `Fatgraph`
        isomorphism is represented.
        """
        return self.isomorphisms(self)


    cpdef list compute_boundary_cycles(self):
        """Return a list of boundary cycles of this `Fatgraph` object.

        Boundary cycles are represented as a cyclic list of 'corners':
        a corner is a triple `(v, i, j)` consisting of a vertex and
        two consecutive indices (in the cyclic order, so, either `j ==
        i+1` or `i` and `j` are the starting and ending indices)::
        
          >>> Fatgraph([Vertex([2,1,0]),Vertex([2,0,1])]).compute_boundary_cycles()
          [BoundaryCycle([(1, 2, 0), (0, 0, 1)]),
           BoundaryCycle([(0, 1, 2), (1, 1, 2)]),
           BoundaryCycle([(0, 2, 0), (1, 0, 1)])]

        This verbose representation allows one to distinguish the
        boundary cycles made from the same set of edges::

          >>> Fatgraph([Vertex([0,1,2,0,1,2])]).compute_boundary_cycles()
          [BoundaryCycle([(0, 2, 3), (0, 4, 5), (0, 0, 1)]),
           BoundaryCycle([(0, 1, 2), (0, 3, 4), (0, 5, 0)])]
        """
        assert self.num_external_edges == 0, \
               "Fatgraph.boundary_cycles: "\
               " cannot compute boundary cycles for" \
               " a graph with nonzero external edges: %s" % self

        cdef int v
        cdef int i
        cdef int j
        cdef int edge
        cdef list corners
        cdef list result
        cdef list triples

        # try cached result first
        if self._cache_bcys is not None:
            return self._cache_bcys

        # Build the collection of "corners" of `graph`,
        # structured just like the set of vertices.
        # By construction, `corners[v][i]` has the the
        # form `(v,i,j)` where `j` is the index following
        # `i` in the cyclic order.
        corners = [ [ (v, i, (i+1)%len(self.vertices[v])) for i in xrange(len(self.vertices[v])) ]
                    for v in xrange(self.num_vertices) ]

        result = list()
        while True:
            # build a list of corners comprising the same boundary
            # cycle: start with one corner, follow the edge starting
            # at the second delimiter of the corner to its other
            # endpoint, and repeat until we come back to the starting
            # point.  
            triples = list()
            (v, i, corner) = first_unused_corner(self.vertices, corners)
            if corner is None:
                break
            start = (v,i)
            while (v,i) != start or len(triples) == 0:
                corners[v][i] = None
                triples.append(corner)
                assert v == corner[0]
                assert i == corner[1]
                j = (corner[2])
                edge = (self.vertices[v][j])
                (v,i) = self.other_end(edge, v, j)
                corner = corners[v][i]
            result.append(BoundaryCycle(triples).set_graph(self))

        self._cache_bcys = result
        return result
        

    cdef bridge(self, int edge1, int side1, int edge2, int side2):
        """Return a new `Fatgraph`, formed by inserting trivalent
        vertices in the middle of edges `edge1` and `edge2` and
        connecting them with a new edge.

          >>> g = Fatgraph([Vertex([0,1,2]), Vertex([0,2,1])])
          >>> g1 = g.bridge(0, 0, 1, 1)
          >>> g1 is g
          False
          >>> g1 == g
          False
          
        Arguments `side1` and `side2` control which side the new edge
        is attached to (valid values are 0 or 1), i.e., which of the
        two inequivalent cyclic orders the new trivalent vertices will
        be given::
        
          >>> g = Fatgraph([Vertex([0,1,2]), Vertex([0,2,1])])
          >>> g1 = g.bridge(0, 0, 1, 0)
          >>> g2 = g.bridge(0, 1, 1, 0)
          >>> g1 == g2
          False

        In more detail: let 0,1,2 be the indices of the edges attached
        to the new vertex in the middle of `edge1`, where 0,1 denote
        the two halves of `edge1`.  If `side1` is `0`, then the new
        trivalent vertex will have the cyclic order [0,1,2]; if
        `side1` is `1`, then 0,1 are swapped and the new trivalent
        vertex gets the cyclic order [1,0,2]::

          >>> g1 == Fatgraph([Vertex([0,1,2]), Vertex([4,2,5]), Vertex([0,4,3]), Vertex([1,5,3])])
          True
          >>> g2 == Fatgraph([Vertex([0,1,2]), Vertex([4,2,5]), Vertex([4,0,3]), Vertex([1,5,3])])
          True

        It is worth noting that this procedure involves 5 edges in
        total, 3 of which need new indices.

          >>> g.num_edges
          3
          >>> g1.num_edges
          6
          
        This function is obviously symmetric: the pairs `edge1, side1`
        and `edge2, side2` can be swapped and the result stays the
        same::

          >>> g3 = g.bridge(1, 0, 0, 0)
          >>> g1 == g3
          True

          >>> g4 = g.bridge(1, 0, 0, 1)
          >>> g2 == g4
          True

        Examples::
        
        1) Bridging different sides of the same edge may yield
        different results::
        
          >>> g.bridge(0, 0, 1, 0)  == Fatgraph([Vertex([0,1,2]), Vertex([4,2,5]), Vertex([0,4,3]), Vertex([1,5,3])])
          True
          
          >>> g.bridge(0, 1, 1, 0) == Fatgraph([Vertex([0,1,2]), Vertex([4,2,5]), Vertex([4,0,3]), Vertex([1,5,3])])
          True

        2) One can connect an edge to itself on different sides::
        
          >>> g.bridge(0, 0, 0, 1) == Fatgraph([Vertex([0,1,2]), Vertex([5,2,1]), Vertex([0,4,3]), Vertex([5,4,3])])
          True

        3) And also with both ends on the same side::
        
          >>> g.bridge(0, 0, 0, 0) == Fatgraph([Vertex([0,1,2]), Vertex([5,2,1]), Vertex([0,4,3]), Vertex([4,5,3])])
          True
          
        """
        assert side1 in [0,1], \
               "Fatgraph.bridge: Invalid value for `side1`: '%s' - should be 0 or 1" % side1
        assert side2 in [0,1], \
               "Fatgraph.bridge: Invalid value for `side2`: '%s' - should be 0 or 1" % side2
        
        cdef int opposite_side1
        cdef int opposite_side2
        cdef int connecting_edge
        cdef int one_half1
        cdef int other_half1
        cdef int one_half2
        cdef int other_half2
        cdef int midpoint1_index
        cdef int midpoint2_index
        cdef int v1a
        cdef int v1b
        cdef int pos1a
        cdef int pos1b
        cdef int v2a
        cdef int v2b
        cdef int pos2a
        cdef int pos2b

        opposite_side1 = (0 if side1 else 1)
        opposite_side2 = (0 if side2 else 1)

        ## assign edge indices
        connecting_edge = (self.num_edges)
        ## break `edge1` in two halves: if `v1a` and `v1b` are the
        ## endpoints of `edge1`, then the "one_half" edge extends from
        ## the `v1a` endpoint of `edge1` to the new vertex
        ## `midpoint1`; the "other_half" edge extends from the
        ## `midpoint1` new vertex to `v1b`.
        one_half1 = (edge1)
        other_half1 = (self.num_edges + 1)
        ## break `edge2` in two halves; if `edge2` is the same edge as
        ## `edge1`, then we are breaking the second half of `edge1` in
        ## two parts.  Otherwise, proceed as above.  In any case, the
        ## "other half" of `edge2` needs a new edge index.
        if edge2 == edge1:
            one_half2 = other_half1
        else:
            one_half2 = edge2
        other_half2 = (self.num_edges + 2)

        ## assign new vertex indices
        midpoint1_index = (self.num_vertices)
        midpoint2_index = (self.num_vertices + 1)

        if side1 == 1:
            midpoint1 = Vertex([other_half1, one_half1, connecting_edge])
        else:
            midpoint1 = Vertex([one_half1, other_half1, connecting_edge])

        if side2 == 1:
            midpoint2 = Vertex([other_half2, one_half2, connecting_edge])
        else:
            midpoint2 = Vertex([one_half2, other_half2, connecting_edge])

        ## two new vertices are added: the mid-points of the connected edges.
        new_vertices = self.vertices + [midpoint1, midpoint2]
        ## the connecting edge has endpoints in the mid-points of
        ## `edge1` and `edge2`, and is *always* in third position.
        new_endpoints_v = self.endpoints_v + [ [midpoint1_index, midpoint2_index] ]
        new_endpoints_i = self.endpoints_i + [ [2,2] ]

        (v1a, v1b) = self.endpoints_v[edge1]
        (pos1a, pos1b) = self.endpoints_i[edge1]
        new_endpoints_v[one_half1] = [v1a, midpoint1_index]
        new_endpoints_i[one_half1] = [pos1a, side1]
        if edge1 != edge2:
            # replace `edge1` with new `other_half1` in the second endpoint
            new_vertices[v1b] = Vertex(new_vertices[v1b][:pos1b]
                                       + [other_half1]
                                       + new_vertices[v1b][pos1b+1:])
            new_endpoints_v.append([midpoint1_index, v1b])  # other_half1
            new_endpoints_i.append([opposite_side1, pos1b]) # other_half1
        else:
            # same edge, "other half" ends at the second endpoint
            new_endpoints_v.append([midpoint1_index, midpoint2_index])
            new_endpoints_i.append([opposite_side1, side2]) 

        # replace `edge2` with new `other_half2` in the second
        # endpoint; again we need to distinguish the special case when
        # `edge1` and `edge2` are the same edge.
        (v2a, v2b) = self.endpoints_v[edge2]
        (pos2a, pos2b) = self.endpoints_i[edge2]
        if edge1 != edge2:
            new_endpoints_v[one_half2] = [v2a, midpoint2_index]
            new_endpoints_i[one_half2] = [pos2a, side2]
        else:
            # `edge1 == edge2`, so `one_half2 == other_half1`
            new_endpoints_v[one_half2] = [midpoint1_index, midpoint2_index]
            new_endpoints_i[one_half2] = [opposite_side1, side2]
        # "other half" of second edge *always* ends at the previous
        # edge endpoint, so replace `edge2` in `v2b`.
        new_vertices[v2b] = Vertex(new_vertices[v2b][:pos2b]
                                   + [other_half2]
                                   + new_vertices[v2b][pos2b+1:])
        new_endpoints_v.append([midpoint2_index, v2b])  # other_half2
        new_endpoints_i.append([opposite_side2, pos2b]) # other_half2

        # build new graph 
        new_edge_numbering = self.edge_numbering + \
                             [other_half1, other_half2, connecting_edge]
        return Fatgraph(new_vertices,
                     endpoints = (new_endpoints_v, new_endpoints_i),
                     num_edges = self.num_edges + 3,
                     num_external_edges = self.num_external_edges,
                     orientation = new_edge_numbering,
                     )
    
    
    cdef int _cmp_orient(self, Fatgraph other, iso):
        cdef Permutation image_edge_numbering
        cdef dict pe
        
        pe = iso[2]
        image_edge_numbering = Permutation([(self.edge_numbering[x],
                                             other.edge_numbering[pe[x]])
                                            for x in xrange(self.num_edges)])
        return image_edge_numbering.sign()
    
    
    #@ocache_weakref
    cdef Fatgraph contract_fg(self, int edgeno):
        """Return new `Fatgraph` obtained by contracting the specified edge.

        Examples::

          >>> Fatgraph([Vertex([2,2,0]), Vertex([0,1,1])]).contract(0)
          Fatgraph([Vertex([1, 1, 0, 0])])
          >>> Fatgraph([Vertex([2,1,0]), Vertex([2,0,1])]).contract(1)
          Fatgraph([Vertex([0, 1, 1, 0])])

        The M_{1,1} trivalent graph yield the same result no matter
        what edge is contracted::

          >>> Fatgraph([Vertex([2,1,0]), Vertex([2,1,0])]).contract(0)
          Fatgraph([Vertex([1, 0, 1, 0])])
          >>> Fatgraph([Vertex([2,1,0]), Vertex([2,1,0])]).contract(1)
          Fatgraph([Vertex([0, 1, 0, 1])])
          >>> Fatgraph([Vertex([2,1,0]), Vertex([2,1,0])]).contract(2)
          Fatgraph([Vertex([1, 0, 1, 0])])
        """
        # check that we are not contracting a loop or an external edge
        assert not self.is_loop(edgeno), \
               "Fatgraph.contract: cannot contract a loop."
        assert (self.endpoints_v[edgeno][0] is not None) \
               and (self.endpoints_v[edgeno][1] is not None), \
               "Fatgraph.contract: cannot contract an external edge."
        assert (edgeno >= 0) and (edgeno < self.num_edges), \
               "Fatgraph.contract: invalid edge number (%d):"\
               " must be in range 0..%d" \
               % (edgeno, self.num_edges)

        cdef int v1
        cdef int v2
        cdef int pos1
        cdef int pos2
        cdef int l1
        cdef int l2
        cdef int x
        cdef int edge
        cdef int cut
        cdef dict renumber_edges
        cdef dict renumber_pos1
        cdef dict renumber_pos2
        cdef dict renumber_vertices
        cdef dict renumber_edge_numbering
        cdef list new_vertices
        cdef list new_endpoints_v
        cdef list new_endpoints_i

        ## Plug the higher-numbered vertex into the lower-numbered one.
        
        # store endpoints of the edge-to-be-contracted
        (v1, v2) = self.endpoints_v[edgeno]
        (pos1, pos2) = self.endpoints_i[edgeno]
        if v1 > v2:
            # swap endpoints so that `v1 < v2`
            v1, v2 = v2, v1
            pos1, pos2 = pos2, pos1

        # save highest-numbered index of vertices to be contracted
        l1 = len(self.vertices[v1]) - 1
        l2 = len(self.vertices[v2]) - 1

        # Build new list of vertices, removing the contracted edge and
        # shifting all indices above:
        #   - edges numbered 0..edgeno-1 are unchanged;
        #   - edges numbered `edgeno+1`.. are renumbered, 
        #     shifting the number down one position;
        #   - edge `edgeno` is kept intact, will be removed by mating
        #     operation (see below).
        renumber_edges = dict([(i+1,i)
                               for i in xrange(edgeno, self.num_edges)])
        # See `ltranslate` in utils.py for how this prescription is
        # encoded in the `renumber_edges` mapping.
        new_vertices = [ Vertex(ltranslate(renumber_edges, v))
                         for v in self.vertices ]

        # Mate endpoints of contracted edge:
        # 1. Rotate endpoints `v1`, `v2` so that the given edge would
        #    appear *last* in `v1` and *first* in `v2` (*Note:* since
        #    the contracted edge has already been deleted and `v1`,
        #    `v2` are *cyclic*, this means that we do the same
        #    operation on `v1` and `v2` alike).
        # 2. Join vertices by concatenating the list of incident
        #    edges;
        # 3. Set new `i1` vertex in place of old first endpoint:
        new_vertices[v1] = Vertex(
            new_vertices[v1][pos1+1:] + new_vertices[v1][:pos1]
            +
            new_vertices[v2][pos2+1:] + new_vertices[v2][:pos2]
            )
        # 4. Remove second endpoint from list of new vertices:
        del new_vertices[v2]

        # vertices with index above `v2` are now shifted down one place
        renumber_vertices = dict([(i+1,i)
                                  for i in xrange(v2, self.num_vertices)])
        # vertex `v2` is mapped to vertex `v1`
        renumber_vertices[v2] = v1
        new_endpoints_v = [ ltranslate(renumber_vertices, ep)
                          for ep in  self.endpoints_v ]
        del new_endpoints_v[edgeno]

        new_endpoints_i = [ ep_i[:] for ep_i in self.endpoints_i ]
        del new_endpoints_i[edgeno]
        # renumber attachment indices, according to the mating of
        # vertices `v1` and `v2`:
        # - on former vertex `v1`:
        #   * indices (pos1+1)..l1 are mapped to 0..(l1-pos1-1)
        #     in the mated vertex;
        #   * index pos1 is deleted;
        #   * indices 0..pos1-1 are mapped to (l1-pos1)..l1-1;
        renumber_pos1 = dict()
        renumber_pos1[pos1] = None
        renumber_pos1.update(dict([(pos1+1+x, x)
                                   for x in xrange(l1-pos1)]))
        renumber_pos1.update(dict([(x, l1-pos1+x)
                                   for x in xrange(pos1)]))
        for edge in self.vertices[v1]:
            if edge == edgeno:
                continue # skip contracted edge
            if v1 == self.endpoints_v[edge][0]:
                new_endpoints_i[renumber_edges.get(edge, edge)][0] = renumber_pos1[self.endpoints_i[edge][0]]
            if v1 == self.endpoints_v[edge][1]:
                new_endpoints_i[renumber_edges.get(edge, edge)][1] = renumber_pos1[self.endpoints_i[edge][1]]
        # - on former vertex `v2`:
        #   * indices (pos2+1)..l2 are mapped to l1..(l1+l2-pos2-1);
        #   * index pos2 is deleted;
        #   * indices 0..pos2-1 are mapped to (l1+l2-pos2)..l1+l2-1:
        renumber_pos2 = dict()
        renumber_pos2[pos2] = None
        renumber_pos2.update(dict([(pos2+1+x, l1+x)
                                   for x in xrange(l2-pos2)]))
        renumber_pos2.update(dict([(x, l1+l2-pos2+x)
                                   for x in xrange(pos2)]))
        for edge in self.vertices[v2]:
            if edge == edgeno:
                continue # skip contracted edge
            if v2 == self.endpoints_v[edge][0]:
                new_endpoints_i[renumber_edges.get(edge, edge)][0] = renumber_pos2[self.endpoints_i[edge][0]]
            if v2 == self.endpoints_v[edge][1]:
                new_endpoints_i[renumber_edges.get(edge, edge)][1] = renumber_pos2[self.endpoints_i[edge][1]]

        ## Orientation of the contracted graph.

        cut = (self.edge_numbering[edgeno])
        renumber_edge_numbering = dict()
        renumber_edge_numbering[cut] = None
        renumber_edge_numbering.update(dict([(x,x) for x in xrange(cut)]))
        renumber_edge_numbering.update(dict([(x+1,x)
                                             for x in xrange(cut,self.num_edges)]))
        new_edge_numbering = [ renumber_edge_numbering[self.edge_numbering[x]]
                         for x in xrange(self.num_edges)
                         if x != edgeno ]
        
        # build new graph 
        return Fatgraph(new_vertices,
                     endpoints = (new_endpoints_v, new_endpoints_i),
                     num_edges = self.num_edges - 1,
                     num_external_edges = self.num_external_edges,
                     orientation = new_edge_numbering,
                     )


    cdef dict edge_orbits(self):
        """Compute orbits of the edges under the action of graph
        automorphism group, and a representative for each orbit.
        
        Returns a dictionary, whose keys are the representatives, and
        whose values are the orbits.  Orbits are represented as Python
        `set` objects.

        Examples::

          >>> Fatgraph([Vertex([0,1,2]), Vertex([0,1,2])]).edge_orbits()
          {0: set([0, 1, 2])}

          >>> Fatgraph([Vertex([1, 0, 2]), Vertex([2, 1, 0])]).edge_orbits()
          {0: set([0, 1, 2])}
          
        """
        # try cached result first
        if self._cache_edge_orbits is not None:
            return self._cache_edge_orbits

        cdef int x
        cdef int y
        cdef dict orbits

        orbits = dict( [(x, set([x])) for x in xrange(self.num_edges)] )
        for a in self.automorphisms():
            edge_permutation = a[2]
            for x in xrange(self.num_edges):
                if x not in orbits:
                    continue
                y = edge_permutation[x]
                if y not in orbits:
                    continue
                # `x` and `y` are in the same orbit, only keep the one
                # with lower abs. value, and remove the other.
                if y > x:
                    orbits[x].update(orbits[y])
                    del orbits[y]
        # check that all elements lie in some orbit
        assert sum([len(set(o)) for o in orbits.itervalues()]) == self.num_edges, \
               "Fatgraph.edge_orbits():" \
               " Computed orbits `%s` do not exhaust edge set `%s`" \
               " [%s.edge_orbits() -> %s]" % (orbits, list(self.edges()), self, orbits)

        self._cache_edge_orbits = orbits
        return orbits


    cdef dict edge_pair_orbits(self):
        """Compute orbits of pairs `(edge1, edge2)` under the action
        of graph automorphism group, and a representative for each
        orbit.
        
        Returns a dictionary, whose keys are the representatives, and
        whose values are the orbits.  Orbits are represented as Python
        `set` objects.

        Examples::

          >>> Fatgraph([Vertex([0,1,2]), Vertex([0,1,2])]).edge_pair_orbits()
          {(0, 1): set([(0, 1), (1, 2), (2, 0)]),
           (0, 0): set([(0, 0), (1, 1), (2, 2)]),
           (0, 2): set([(1, 0), (0, 2), (2, 1)])}
          
        """
        # try cached result first
        if self._cache_edge_pair_orbits is not None:
            return self._cache_edge_pair_orbits
        
        cdef int x
        cdef int y
        cdef list edge_pairs
        cdef dict orbits
        
        edge_pairs = [ (x,y) 
                       for x in xrange(self.num_edges)
                       for y in xrange(self.num_edges) ]
        orbits = dict( [(p, set([p])) for p in edge_pairs] )
        for a in self.automorphisms():
            edge_permutation = a[2]
            for p in edge_pairs:
                if p not in orbits:
                    continue
                q = (edge_permutation[p[0]], edge_permutation[p[1]])
                if q not in orbits:
                    continue
                # `p` and `q` are in the same orbit, only keep the one
                # with lower abs. value, and remove the other.
                if p < q:
                    orbits[p].update(orbits[q])
                    del orbits[q]
        # check that all elements lie in some orbit
        assert sum([len(set(o)) for o in orbits.itervalues()]) == len(edge_pairs), \
               "Fatgraph.edge_pair_orbits():" \
               " Computed orbits `%s` do not exhaust edge pairs set `%s`" \
               " [%s.edge_pair_orbits() -> %s]" % (orbits, edge_pairs, self, orbits)

        self._cache_edge_pair_orbits  = orbits
        return orbits


    cdef endpoints(self, int edgeno):
        """Return the endpoints of `edge`, as a pair of `(v, pos)`
        where `v` is the endpoint vertex index, and `pos` is the
        attachment index of `edge` into the `Vertex` object
        `self.vertices[v]`.

        The pair `((v1, pos1), (v2, pos2))` is ordered such that `v1 < v2`.
        """
        cdef int v1
        cdef int v2
        cdef int pos1
        cdef int pos2
        
        v1, v2 = self.endpoints_v[edgeno]
        pos1, pos2 = self.endpoints_i[edgeno]
        if v1 > v2:
            # swap endpoints so that `v1 < v2`
            v1, v2 = v2, v1
            pos1, pos2 = pos2, pos1
        return ((v1, pos1), (v2, pos2))


    cdef Fatgraph hangcircle(self, int edge, int side):
        """Return a new `Fatgraph`, formed by attaching a circle with
        a new edge to a new trivalent vertex in the middle of `edge`.

          >>> g = Fatgraph([Vertex([0,1,2]), Vertex([0,2,1])])
          >>> g1 = g.hangcircle(0, 0)
          >>> g1 is g
          False
          
        Argument `side` controls which side of `edge` the circle is
        hung to (valid values are 0 or 1), i.e., which of the two
        inequivalent cyclic orders the new trivalent vertices will be
        given::
        
          >>> g = Fatgraph([Vertex([0,1,2]), Vertex([0,2,1])])
          >>> g1 = g.hangcircle(0, 0)
          >>> g1 == Fatgraph([Vertex([0,1,2]), Vertex([3,2,1]), Vertex([0,3,4]), Vertex([5,5,4])])
          True
          >>> g2 = g.hangcircle(0, 1)
          >>> g2 == Fatgraph([Vertex([0,1,2]), Vertex([3,2,1]), Vertex([3,0,4]), Vertex([5,5,4])])
          True

        It is worth noting that the new graph will have 3 edges more
        than the original one::

          >>> g1.num_edges == g.num_edges + 3
          True
          
        """
        assert side in [0,1], \
               "Fatgraph.hangcircle: Invalid value for `side`: '%s' - should be 0 or 1" % side
        
        cdef int opposite_side
        cdef int one_half
        cdef int other_half
        cdef int connecting_edge
        cdef int circling_edge
        cdef int midpoint_index
        cdef int T_index
        cdef int v1
        cdef int v2
        cdef int pos1
        cdef int pos2
        cdef list new_vertices
        cdef list new_endpoints_v
        cdef list new_endpoints_i
        cdef list new_edge_numbering

        opposite_side = (0 if side else 1)

        ## assign edge indices
        
        ## break `edge` in two halves: if `v1` and `v2` are the
        ## endpoints of `edge`, then the "one_half" edge extends from
        ## the `v1` endpoint of `edge1` to the new vertex
        ## `midpoint`; the "other_half" edge extends from the
        ## `midpoint` new vertex to `v2`.
        one_half = (edge)
        other_half = (self.num_edges)

        connecting_edge = (self.num_edges + 1)
        circling_edge = (self.num_edges + 2)
        
        ## assign new indices to new vertices
        midpoint_index = (self.num_vertices)
        T_index = (self.num_vertices + 1)

        ## two new vertices are added: the mid-point of `edge`, and
        ## the vertex `T` lying on the circle.
        if side == 1:
            midpoint = Vertex([other_half, one_half, connecting_edge])
        else:
            midpoint = Vertex([one_half, other_half, connecting_edge])
        T = Vertex([circling_edge, circling_edge, connecting_edge])
        new_vertices = self.vertices + [midpoint, T]

        ## new edge endpoints:
        ## - start with a copy of the original ednpoints:
        new_endpoints_v = self.endpoints_v[:]
        new_endpoints_i = self.endpoints_i[:]

        ## - replace `edge` with new `other_half` in the second endpoint:
        (v1, v2) = self.endpoints_v[edge]
        (pos1, pos2) = self.endpoints_i[edge]
        new_endpoints_v[one_half] = [v1, midpoint_index]
        new_endpoints_i[one_half] = [pos1, side]
        new_vertices[v2] = Vertex(new_vertices[v2][:pos2]
                                  + [other_half]
                                  + new_vertices[v2][pos2+1:])
        new_endpoints_v.append([midpoint_index, v2])  # other_half1
        new_endpoints_i.append([opposite_side, pos2]) # other_half1

        ## - the connecting edge has endpoints in the mid-point of
        ## `edge` and in `T`, and is *always* in third position:
        new_endpoints_v.append([midpoint_index, T_index])
        new_endpoints_i.append([2,2])

        ## - the circling edge is a loop with vertex `T`
        new_endpoints_v.append([T_index, T_index])
        new_endpoints_i.append([0,1])
        
        # finally, build new graph 
        new_edge_numbering = self.edge_numbering + \
                             [other_half, connecting_edge, circling_edge]
        return Fatgraph(new_vertices,
                     endpoints = (new_endpoints_v, new_endpoints_i),
                     num_edges = self.num_edges + 3,
                     num_external_edges = self.num_external_edges,
                     orientation = new_edge_numbering,
                     )


    cpdef bint is_loop(self, int edge):
        """Return `True` if `edge` is a loop (i.e., the two endpoint coincide).
        """
        return self.endpoints_v[edge][0] == self.endpoints_v[edge][1]
        

    cpdef bint is_orientation_reversing(self, automorphism):
        """Return `True` if `automorphism` reverses orientation of
        this `Fatgraph` instance."""
        return (-1 == Fatgraph._cmp_orient(self, self, automorphism))


    cpdef bint is_oriented(self):
        """Return `True` if `Fatgraph` is orientable.

        A `Fatgraph` is orientable iff it has no orientation-reversing
        automorphism.

        Enumerate all automorphisms, end exit with `False` result as
        soon as one orientation-reversing one is found.

        Examples::

          >>> Fatgraph([Vertex([1,0,1,0])]).is_oriented()
          False

          >>> Fatgraph([Vertex([2, 0, 1]), Vertex([2, 0, 1])]).is_oriented()
          True
          
          >>> Fatgraph([Vertex([2, 1, 0]), Vertex([2, 0, 1])]).is_oriented()
          False
          
          >>> Fatgraph([Vertex([2, 1, 1]), Vertex([2, 0, 0])]).is_oriented()
          False

          >>> Fatgraph([Vertex([3, 2, 2, 0, 1]), Vertex([3, 1, 0])], \
                    numbering=[(0, CiclicList((2,))),  \
                               (1, CiclicList((0, 1))),  \
                               (2, CiclicList((3, 1))),  \
                               (3, CiclicList((0, 3, 2))) ]) \
                               .is_oriented()
          True
          >>> Fatgraph([Vertex([2, 3, 1]), Vertex([2, 1, 3, 0, 0])], \
                       numbering=[(0, CiclicList((0,))),  \
                                  (2, CiclicList((1, 3))),  \
                                  (3, CiclicList((3, 0, 2))),  \
                                  (1, CiclicList((2, 1))) ]) \
                               .is_oriented()
          True
        """
        ## Try to find an orientation-reversing automorphism the hard way
        for a in self.automorphisms():
            if self.is_orientation_reversing(a):
                return False
        # no orientation reversing automorphism found
        return True


    cpdef list isomorphisms(g1, g2):
        """Return list of `Fatgraph` isomorphisms from `g1` to `g2`.

        An isomorphism is represented by a tuple `(pv, rot, pe)` where:

          - `pv` is a permutation of ther vertices: the `i`-th vertex
            of `g1` is sent to the `pv[i]`-th vertex of `g2`, rotated
            by `rot[i]` places leftwards;

          - `pe` is a permutation of the edges: edge `i` in `g1` is
            mapped to edge `pe[i]` in `g2`.

        This method can iterate over the automorphism group of a
        graph::

          >>> g1 = Fatgraph([Vertex([2, 1, 1]), Vertex([2, 0, 0])])
          >>> for f in g1.isomorphisms(g1): print f
          ({0: 0, 1: 1}, [0, 0], {0: 0, 1: 1, 2: 2})
          ({0: 1, 1: 0}, [0, 0], {0: 1, 1: 0, 2: 2})

        Or it can find the isomorphisms between two given graphs::

          >>> g2 = Fatgraph([Vertex([2, 2, 0]), Vertex([1, 1, 0])])
          >>> for f in g1.isomorphisms(g2): print f
          ({0: 0, 1: 1}, [2, 2], {0: 1, 1: 2, 2: 0})
          ({0: 1, 1: 0}, [2, 2], {0: 2, 1: 1, 2: 0})

        If there are no isomorphisms connecting the two graphs, then no
        item is returned by the iterator::

          >>> g3 = Fatgraph([Vertex([2, 1, 0]), Vertex([2, 0, 1])])
          >>> g1.isomorphisms(g3)
          []
        """
        cdef int val
        cdef int src0
        cdef int dst0
        cdef int rot0
        cdef int x
        cdef dict vs1
        cdef dict vs2
        cdef Permutation pv
        cdef list rot
        cdef Permutation pe
        cdef list nexts
        cdef list neighborhood

        cdef object key = g2.cache_id()
        try:
            return g1._cache_isomorphisms[key]
        except KeyError:
            pass # continue with building isomorphisms list
        
        # As this procedure is quite complex, we break it into a
        # number of auxiliary functions.
        
        # if graphs differ in vertex valences, no isomorphisms
        vs1 = g1.valence_spectrum()
        vs2 = g2.valence_spectrum()
        if not set(vs1.keys()) == set(vs2.keys()):
            return list()
        # if graphs have unequal vertex distribution by valence, no isomorphisms
        for val in g1.vertex_valences():
            if len(vs1[val]) != len(vs2[val]):
                return list()

        if g1._cache_starting_vertices is None:
            g1._cache_starting_vertices = starting_vertices(g1)
        (val, vs) = g1._cache_starting_vertices
        src0 = (vs[0])
        v = g1.vertices[src0]
        result = list()
        for dst0 in admissible_vertex_mappings(v, g2, vs2[val]):
            for rot0 in xrange(val):
                try:                    
                    # pass 0: init new (pv, rot, pe) triple
                    pv = Permutation()
                    rot = [ None for x in xrange(g1.num_vertices) ]
                    pe = Permutation()

                    # pass 1: map `v` to `v2` and build map
                    # of neighboring vertices for next pass
                    pv[src0] = dst0
                    rot[src0] = rot0
                    if not pe.extend(v, g2.vertices[dst0][rot0:rot0+val]):
                        continue # to next `rot0`

                    # pass 2: extend map to neighboring vertices
                    m = (pv, rot, pe)
                    nexts = neighbors(m, g1, src0, g2, dst0)
                    while len(m[0]) < g1.num_vertices:
                        neighborhood = list()
                        for (i1, i2, r) in nexts:
                            m = extend_map(m, g1, i1, r, g2, i2)
                            neighborhood += neighbors(m, g1, i1, g2, i2)
                        nexts = neighborhood

                # extension failed in the above block, continue with next candidate
                except CannotExtendMap:
                    continue # to next `rot0`

                # sanity checks
                assert len(m[0]) == g1.num_vertices
                assert len(m[2]) == g1.num_edges
                assert set(m[0].keys()) == set(range(g1.num_vertices))
                assert set(m[0].values()) == set(range(g2.num_vertices))
                assert set(m[2].keys()) == set(range(g1.num_edges))
                assert set(m[2].values()) == set(range(g2.num_edges))

                # finally
                result.append(m)

        # remember result
        g1._cache_isomorphisms[key] = result

        # at last
        return result


    cpdef int num_automorphisms(self):
        """Return the cardinality of the automorphism group of this
        `Fatgraph` object.

        Examples::

          >>> Fatgraph([Vertex([0,1,2]), Vertex([0,2,1])]).num_automorphisms()
          6
          >>> Fatgraph([Vertex([0,1,1]), Vertex([0,2,2])]).num_automorphisms()
          2
        """
        return len(self.automorphisms())
    
        
    # auxiliary function
    cdef other_end(graph, int edge, vertex, int attachment):
        """Return the other endpoint of `edge` as a pair `(v, i)`.
        """
        ends = zip(graph.endpoints_v[edge], graph.endpoints_i[edge])
        if ends[0] == (vertex, attachment):
            return ends[1]
        else:
            return ends[0]


    cpdef int projection(self, other):
        """Return the component of the projection of `self` on the
        basis vector `other`.  This can be either 0 (if `self` and
        `other` are not isomorphic), or +1/-1 depending on comparison
        of the orientation of `self` with the pull-back orientation on
        `other`.

        If the two graphs are *not* isomorphic, then the result is 0::

          >>> g1 = Fatgraph([Vertex([0,1,2,0,1,2])])
          >>> g2 = Fatgraph([Vertex([0,1,2]), Vertex([0,1,2])])
          >>> Fatgraph.projection(g1, g2)
          0

        Any graph obviously projects onto itself with coefficient `1`::

          >>> Fatgraph.projection(g1, g1)
          1
          >>> Fatgraph.projection(g2, g2)
          1

        And similarly does any graph isomorphic to a given graph::
        
          >>> g3 = Fatgraph(g2)
          >>> Fatgraph.projection(g2, g3)
          1

        Flipping the orientation on an edge reverses the coefficient
        sign::

          >>> g4 = Fatgraph(g2)
          >>> # make a copy of `g4.edge_numbering`, since it's shared with `g2`
          >>> g4 .edge_numbering = copy(g4.edge_numbering)
          >>> g4.edge_numbering[0], g4.edge_numbering[1] = \
              g4.edge_numbering[1], g4.edge_numbering[0]
          >>> Fatgraph.projection(g2, g4)
          -1
          
        """
        assert isinstance(other, Fatgraph), \
               "Fatgraph.projection:" \
               " called with non-Fatgraph argument `other`: %s" % other
        assert self.is_oriented(), \
               "Fatgraph.projection: cannot project non-orientable graph: %s" \
               % self
        assert other.is_oriented(), \
               "Fatgraph.projection: cannot project non-orientable graph: %s" \
               % other
        iso = self.isomorphisms(other)
        if len(iso) > 0:
            # since both graphs are orientable, it does not matter
            # which isomorphism we pick: if `f` and `f'` are two
            # isomorphisms, then `f' = a' f a` where `a` and `a'` are
            # automorphisms of `g1` and `g2`, hence `f` and `f'` have
            # the same sign.
            return Fatgraph._cmp_orient(self, other, iso[0])
        else:
            # list of morphisms is empty, graphs are not equal.
            return 0

    cpdef dict valence_spectrum(self):
        """Return a dictionary mapping valences into vertex indices.

        Examples::

           >>> Fatgraph([Vertex([1,1,0,0])]).valence_spectrum()
           {4: [0]}

           >>> Fatgraph([Vertex([1,1,0]), Vertex([2,2,0])]).valence_spectrum()
           {3: [0, 1]}

           >>> Fatgraph([Vertex([3, 1, 0, 1]), \
                      Vertex([4, 4, 0]), Vertex([3, 2, 2])]).valence_spectrum()
           {3: [1, 2], 4: [0]}
        """
        # try cached result first
        if self._cache_valence_spectrum is not None:
            return self._cache_valence_spectrum
        
        cdef int l
        cdef dict result
        
        result = dict()
        for (index, vertex) in enumerate(self.vertices):
            l = (len(vertex))
            if l in result:
                result[l].append(index)
            else:
                result[l] = [index]
        # consistency checks
        assert set(result.keys()) == set(self.vertex_valences()), \
               "Fatgraph.valence_spectrum:" \
               "Computed valence spectrum `%s` does not exhaust all " \
               " vertex valences %s" \
               % (result, self.vertex_valences())
        assert set(lconcat(result.values())) \
               == set(range(self.num_vertices)), \
               "Fatgraph.valence_spectrum:" \
               "Computed valence spectrum `%s` does not exhaust all " \
               " %d vertex indices" % (result, self.num_vertices)

        self._cache_valence_spectrum = result
        return result

    #@ocache0
    cpdef frozenset vertex_valences(self):
        return frozenset([len(v) for v in self.vertices])

    #@ocache0
    cpdef dict vertex_valence_distribution(self):
        spec = self.valence_spectrum()
        return dict([(v, len(spec[v]))
                     for v in spec.iterkeys()])
    

# auxiliary function in `MgnTrivalentGraphsRecursiveGenerator`
cdef inline bint _admissible(Fatgraph G, int g, int n):
    """Return `True` if `Fatgraph` `G` has genus `g` and number of
    boundary cycles `n`.

    Examples::

      >>> G1 = Fatgraph([Vertex([0,1,2]), Vertex([0,2,1])])
      >>> _admissible(G1, 0, 3)
      True
      >>> _admissible(G1, 1, 1)
      False

      >>> G2 = Fatgraph([Vertex([0,1,2]), Vertex([0,1,2])])
      >>> _admissible(G2, 0, 3)
      False
      >>> _admissible(G2, 1, 1)
      True
    """
    return (G.genus, G.num_boundary_cycles) == (g,n)

cpdef list MgnTrivalentGraphsRecursiveGenerator(int g, int n):
    """Return a list of all connected trivalent fatgraphs having the
    prescribed genus `g` and number of boundary cycles `n`.
    
    Examples::

      >>> for g in MgnTrivalentGraphsRecursiveGenerator(0,3): print g
      Fatgraph([Vertex([1, 2, 1]), Vertex([2, 0, 0])]) 
      Fatgraph([Vertex([1, 0, 2]), Vertex([2, 0, 1])])

      >>> for g in MgnTrivalentGraphsRecursiveGenerator(1,1): print g
      Fatgraph([Vertex([1, 0, 2]), Vertex([2, 1, 0])])

    """
    cdef Fatgraph G
    cdef Fatgraph G_
    cdef int x
    cdef int y
    cdef list unique
    cdef object checkpoint
    cdef object error

    # avoid infinite recursion in later statements
    if n==0 or (g,n)<(0,3):
        return list()

    # sanity check
    assert n > 0, \
           "MgnTrivalentGraphsRecursiveGenerator: " \
           " number of boundary cycles `n` must be positive,"\
           " but got `%s` instead" % n

    logging.debug("Starting MgnTrivalentGraphsRecursiveGenerator(%d,%d) ..." % (g,n))

    ## M_{0,3} - induction base
    if (g,n) == (0,3):
        logging.debug("  MgnTrivalentGraphsRecursiveGenerator(0,3) done.")
        return list([Fatgraph([Vertex([1, 2, 1]), Vertex([2, 0, 0])]),
                     Fatgraph([Vertex([1, 0, 2]), Vertex([2, 0, 1])])])

    ## M_{1,1} - induction base
    elif (g,n) == (1,1):
        logging.debug("  MgnTrivalentGraphsRecursiveGenerator(1,1) done.")
        return list([Fatgraph([Vertex([1, 0, 2]), Vertex([2, 1, 0])])])

    ## General case
    else:
        unique = None
        checkpoint = (runtime.options.checkpoint_dir + '/'
                      + "MgnTrivalentGraphsRecursiveGenerator%d,%d.txt"
                      % (g, n))
        if runtime.options.restart:
            try:
                unique = load(checkpoint)
                logging.debug("  Loaded %d trivalent graphs from file '%s'",
                              len(unique), checkpoint)
            except Exception, error:
                logging.debug("  Could not load saved state from file '%s': %s",
                              checkpoint, error.message)
                # ignore error and proceed

        if unique is None:
            unique = list()

            logging.debug("  MgnTrivalentGraphsRecursiveGenerator(%d,%d): "
                          "pass 1: hang a circle to all edges of graphs in M_{%d,%d} ..." % (g,n, g,n-1))
            for G in MgnTrivalentGraphsRecursiveGenerator(g,n-1):
                for x in G.edge_orbits():
                    G_ = G.hangcircle(x,0)
                    if _admissible(G_, g, n) and (G_ not in unique):
                        unique.append(G_)
                    G_ = G.hangcircle(x,1)
                    if _admissible(G_, g, n) and (G_ not in unique):
                        unique.append(G_)

            logging.debug("  MgnTrivalentGraphsRecursiveGenerator(%d,%d): "
                          "pass 2: bridge all edges of a single graph in M_{%d,%d} ..." % (g,n, g,n-1))
            for G in MgnTrivalentGraphsRecursiveGenerator(g,n-1):
                for (x,y) in G.edge_pair_orbits():
                        G_ = G.bridge(x,0, y,0)
                        if _admissible(G_, g, n) and (G_ not in unique):
                            unique.append(G_)
                        G_ = G.bridge(x,0, y,1)
                        if _admissible(G_, g, n) and (G_ not in unique):
                            unique.append(G_)
                        G_ = G.bridge(x,1, y,0)
                        if _admissible(G_, g, n) and (G_ not in unique):
                            unique.append(G_)
                        G_ = G.bridge(x,1, y,1)
                        if _admissible(G_, g, n) and (G_ not in unique):
                            unique.append(G_)

            logging.debug("  MgnTrivalentGraphsRecursiveGenerator(%d,%d): "
                          "pass 3: bridge all edges of a single graph in M_{%d,%d} ..." % (g,n, g-1,n+1))
            for G in MgnTrivalentGraphsRecursiveGenerator(g-1,n+1):
                for (x,y) in G.edge_pair_orbits():
                        G_ = G.bridge(x,0, y,0)
                        if _admissible(G_, g, n) and (G_ not in unique):
                            unique.append(G_)
                        G_ = G.bridge(x,0, y,1)
                        if _admissible(G_, g, n) and (G_ not in unique):
                            unique.append(G_)
                        G_ = G.bridge(x,1, y,0)
                        if _admissible(G_, g, n) and (G_ not in unique):
                            unique.append(G_)
                        G_ = G.bridge(x,1, y,1)
                        if _admissible(G_, g, n) and (G_ not in unique):
                            unique.append(G_)

            # clear cached isomorphisms, they will not be needed in the sequel
            for G in unique:
                G._cache_isomorphisms.clear()

            # save to checkpoint file
            save(unique, checkpoint)
            
            logging.debug("  MgnTrivalentGraphsRecursiveGenerator(%d,%d) done: %d unique graphs." % (g,n, len(unique)))
            
        return unique



cdef class MgnGraphsIterator(BufferingIterator):
    """Iterate over all connected fatgraphs having the
    prescribed genus `g` and number of boundary cycles `n`.
    
    Examples::

      >>> for g in MgnGraphsIterator(0,3): print g
      Fatgraph([Vertex([1, 2, 1]), Vertex([2, 0, 0])]) 
      Fatgraph([Vertex([1, 0, 2]), Vertex([2, 0, 1])])
      Fatgraph([Vertex([1, 1, 0, 0])])

      >>> for g in MgnGraphsIterator(1,1): print g
      Fatgraph([Vertex([1, 0, 2]), Vertex([2, 1, 0])])
      Fatgraph([Vertex([1, 0, 1, 0])])

    """

    def __init__(self, g, n):
        assert n > 0, \
               "MgnGraphsIterator: " \
               " number of boundary cycles `n` must be positive,"\
               " but got `%s` instead" % n
        assert (g > 0) or (g == 0 and n >= 3), \
               "MgnGraphsIterator: " \
               " Invalid (g,n) pair (%d,%d): "\
               " need either g>0 or g==0 and n>2" \
               % (g,n)

        #: Prescribed genus of returned graphs
        self.g = g

        #: Prescribed number of boundary components
        self.n = n
        
        # Gather all 3-valent graphs.
        cdef list trivalent
        trivalent = list(MgnTrivalentGraphsRecursiveGenerator(g,n))

        #: Fatgraphs to be contracted at next `.refill()` invocation
        self._batch = trivalent

        #: Graphs returned by next `.refill()` call will have this
        #  number of vertices.
        self._num_vertices = 4*g + 2*n - 5
        
        # initialize superclass with list of trivalent graphs
        BufferingIterator.__init__(self, trivalent)


    cdef list refill(self):
        if self._num_vertices == 0:
            raise StopIteration

        cdef int discarded
        cdef int edge
        cdef Fatgraph graph
        cdef Fatgraph dg
        cdef list next_batch
        cdef object checkpoint
        cdef object error

        checkpoint = (runtime.options.checkpoint_dir + '/' + "M%d,%d-MgnGraphsIterator%d.txt"
                      % (runtime.g, runtime.n, self._num_vertices))

        next_batch = None
        if runtime.options.restart:
            try:
                next_batch = load(checkpoint)
                logging.info("  Loaded graphs with %d vertices from file '%s'",
                             self._num_vertices, checkpoint)
            except Exception, error:
                logging.debug("  Could not retrieve state from file '%s': %s",
                              checkpoint, error.message)
                # ignore error, continue generating graphs

        if next_batch is None:
            logging.debug("Generating graphs with %d vertices ...",
                         self._num_vertices)
            next_batch = list()
            discarded = (0)
            for graph in self._batch:
                # contract all edges
                for edge in graph.edge_orbits():
                    if not graph.is_loop(edge):
                        dg = graph.contract_fg(edge)
                        if dg not in next_batch:
                            # put graph back into next batch for processing
                            next_batch.append(dg)
                        else:
                            dg.release()
                            discarded += 1
            logging.info("  Found %d distinct unique fatgraphs with %d vertices, discarded %d duplicates.",
                         len(self._batch), self._num_vertices, discarded)

            # clear cached isomorphisms, they will not be needed
            for G in next_batch:
                G._cache_isomorphisms.clear()

            # save generated graphs to checkpoint file
            save(next_batch, checkpoint)

        self._batch = next_batch
        self._num_vertices -= 1
        return next_batch

#MgnGraphsIterator = persist.PersistedIterator(_MgnGraphsIterator)



## main: run tests

if "__main__" == __name__:
    import doctest
    doctest.testmod(name='rg',
                    optionflags=doctest.NORMALIZE_WHITESPACE)
